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Abstract 

A novel integral equations approach is applied for studying ion pairing in the restricted prim- 
itive model (RPM) electrolyte, i. e., the three point extension (TPE) to the Ornstein-Zernike 
integral equations. In the TPE approach, the three-particle correlation functions (ri,r2,r3) 
are obtained. The TPE results are compared to molecular dynamics (MD) simulations and other 
theories. Good agreement between TPE and MD is observed for a wide range of parameters, par- 
ticularly where standard integral equations theories fail, i. e., low salt concentration and high ionic 
valence. Our results support the formation of ion pairs and aligned ion complexes. 
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I. INTRODUCTION 



The restricted primitive model electrolyte (RPM) has been widely studied by means of 



integral equations 



flflflUflflflS 



and Monte Carlo (MC) simulations 
lid phase 

IS Eli 



101. All the 



approaches describe well the RPM in a wide regime of the fluid phase diagram. Nonetheless, 
they all fail in the dilute regime of a multivalent electrolyte J, Q, which can be relevant 



for the study of phase transitions in ionic fluids. 

Such phase transition have been predicted as early as 1962, for ionic mixtures 



later for polyelectrolyte solutions 



121, and 



131 ] - Experiments for the gas-liquid phase transition of 



molten salts have been made in the past 



14j . Among the first computer simulations of the 



RPM, where this transition is reported are those of Vorontsov-Vel'yaminov et al. jla, |R 
Recently, there has been a renewed interest in ionic phase transitions: Computer simulations 

18 . 3, and for variations of this model, where unsymmetrical ionic 



studies for the RPM 



17 



charge and size is considered 



20 



21 



22 



23, 



MUm, have given insight into the nature of the 



phase transition and the molecular mechanisms behind these transitions. Experiments of a 



26 



27 



28, 



29 



Phase transitions 



liquid-liquid phase transitions have also been reported 
of the RPM can be identified either with molten salts gas-liquid transition or with the two 
liquid transition, since in terms of dimensionless parameters the RPM does not distinguish 
between these two sistems |3lj |. 

In the past, it has been proposed a powerful approach to systematically incorporate 



Hveriui approacn to 5 



This method is known as 



correlations into any given liquid theory 
three point extension (TPE) to integral equations. By construction, TPE explicitly provides 
valuable information of the three-particle correlations. In consequence, the resulting pair 
distribution function includes virtually infinitely more correlations and, hence, a better 
system description is expected. Although in the past TPE has not been applied to bulk 
fluids, our presumption is sustained by previous TPE calculation for inhomogeneous fluids 
where better agreement with computer simulations S. 3ol . I^J were reported than in the 
case of standard j^l] integral equation theories [3, [3]. 

In this paper we apply TPE to the RPM and compare with our molecular dynamics 
(MD) simulations. Based on ion yairinq association, phase transitions in ionic fluids have 



been reported by computer simulations 



24j . The results presented here, (both TPE and our 



MD) for this region of the phase diagram, support this ion paring association mechanisms. 
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Moreover, TPE, based on the agreement with our MD results, provides a reliable theory to 
study ionic fluids in the important phase diagram region of low ionic concentration and high 
coulombic coupling. 

In spite of important theoretical efforts made in the past, a proper description of the full 
RPM electrolyte phase diagram is still required. Previous approaches to study triplet corre- 



lations have been developed by Kjellander et al. |42j and Plischke and Henderson [43, 44]. 



In their study, they considered a fluid next to a plate and they computed the inhomogeneous 



two-partic 



e distribution function. More relevant for the present study, however, is that of 



Attard 45j, |46j , who calculated the two particle inhomogeneous distribution function (using 
the Percus-Yevick closure) for a hard sphere fluid next to a hard sphere particle. In his 
approach, he finds an excellent agreement with MC data. To the best of our knowledge, no 
triplet correlation function has been explicitly calculated for the RPM electrolyte. 

In a study of the critical behavior of the RPM electrolyte at the level of the Debye- 



Hiickel theory, Levin and Fisher 



47 



have included triplet correlations by imposing the 
presence of ionic pairs, such a consideration reveals an Ising critical behavior. The ions pairs 
idea first proposed by Bierr um 14911 has been considerably extende d by Levin and Fisher, 

and Blum and co-workers 



51 



52, 5c 



d by Levin and 
3H 57, Q- 



While 



Stell and co-workers 

ion pairing clearly seems to be the molecular mechanism ruling the ionic solutions phase 



transitions 



18. 



21 



24. 



28, 



its physical bases remain unexplained 



241. On the other 

nn 



hand, although in the past s 0m e e X peri m ents supported a dassica, critical behavior t3t tad, 
others the Ising universality class [61] and Singh and Pitzer [62] suggested a crossover from 
classical to I sing behavior, later experimental results, however, seem to agree in a crossover 



behavior 



3, 



30 



63 



64 



651 ] . Therefore, while ionic fluids asymptotic critical behavior appear 



to exhibit ultimately Ising-like critical behavior, the question of why do some ionic fluids 
appear to display classical behavior j^, remains unanswered 6^- A shortcoming of the ion 
pairing theories is that ion pairing is imposed and hence they provide an ad hoc molecular 
mechanism. Perhaps a better understanding of the molecular mechanisms behind phase 
transitions could be captured by a formal many body theory, such as TPE, where three 
particle correlations are calculated explicitly, and no ion pairing is imposed. 

In this work, by using the TPE to integral equations approach, we obtain a better de- 
scription of the RPM electrolyte: In particular for the strongly coupled region. We also 
analyze the formation of ion complexes. The structure of the article is set out as follows. In 
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Sec. |n]we present the TPE formalism. Section HTT1 is devoted to the computational details of 
the MD simulation. In Sec. IIVI we present our results for two typical (divalent) electrolyte 
concentrations. The obtained three particle distribution function, (r 1; r 2 , r 3 ), with TPE 
and MD simulation are compared and analyzed in terms of ion asociation. We also compare 
the mean force between two particles obtained with TPE, conventional HNC/MSA and MD. 
Finally, Sec. IS contains concluding remarks. 



II. THEORY 



A. Three point extension to integral equation theories 

The pair correlation function, g(ri 2 = i"i — r2), of a one-component fluid with its compo- 
nents interacting through the pair potential ii(r 12 ), is related to the potential of mean force 
w(ri 2 ) (between two particles located at ri and r 2 ) by 

g(r 12 ) = exp{-/3w(r 12 )} . (1) 

If g (i"i 2 ) is expanded in powers of the bulk concentration, the n-th order coefficient is a 
sum of integrals of products of the Mayer function /(r 12 ) = exp{— /3w(r 12 )} — 1. Such an 
integral of a product of Mayer functions can be conveniently represented by Mayer diagrams 
js?! IqtI loT^ . The diagrams of the first and second order coefficients are given in the left hand 
side of Fig. ^ There is still not an exact theory to compute g(r 12 ), and all the available 
theories ignore several classes of topologically different diagrams. We will come back to 
this point below when we discuss the direct correlation function and the Ornstein-Zernike 
equation. 

In a mult i- component fluid, the total correlation function, hij (ri 2 ) = g%j (1*12) — 1, between 
two particles of species i and j located at ri and r 2 , respectively, is related to the direct cor- 
relation function, Cy (i"i 2 ), through the Ornstein-Zernike equation which for a ^-component 
fluid is given by 

hij (r 12 ) = Cij (r 12 ) + (1*23) c m j (i"i3) dr 3 , (2) 

771=1 

where p m is the concentration of species m. Several closures between (r 12 ) and (r 12 ) 
have been proposed. For instance, 



Cij (ri 2 ) = -(3iiij (ri 2 ) + % (r i2 ) - ln&j (ri 2 ) , 



(3) 



Cij (ri 2 ) = (r i2 ) , and (4) 

Oi 2 ) = /ij (r 12 ) g tj (r 12 ) exp {/fay (r 12 )} . (5) 

Equations (j3J), (JU) and (0) are known as the hypernetted chain equation (HNC), the mean 
spherical approximation (MSA) and the Percus-Yevick (PY) equation, respectively. In the 
hypernetted chain theory, the bridge diagrams are ignored whereas in the Percus-Yevick 
approximation both the bridge and product diagrams are neglected [66, 68]. The first and 
second order Mayer graphs of the HNC and PY theories are also given in Fig. ^ 



Let us now propose 



, |37( that in a fluid of k -species there is an additional dumbbell 
species at infinite dilution made up of two particles (of species f3 and 7) at fixed relative 
position t = r 12 (see Fig. EJ). By defining the dumbbell species as a, we now have a (k + 1) 
-component fluid. For p a — > 0, the total correlation function between the particle of species 
a and the fluid particle of species j reads 

Kj (r 3 ) = c aj (r 3 ) + Y, Pm Km (r 4 ) c mj (r 34 ) dr 4 , (6) 
m=l 

where c m j (r^) is the direct correlation function between particles of species m and j both 
different from a. In order to obtain c m j{x^i), the fc-component Ornstein-Zernike equation 
[Eq. (J2J)] has to be used. Different integral equation theories can be obtained depending 
on the closure relations used for c a j (r 3 ) and c m j (r 34 ) in Eq. ©. For instance, TPE- 
HNC/MSA is obtained if MSA [Eq. ©] is used for c mj (r 34 ) and HNC [Eq.©] for c ai (r 3 ). 

In this formalism, the distribution function, g ai (r 3 ), of the i species around the a 
species can be interpreted as a conditional three-particle distribution function denoted by 
9pii( r 3'i^) = ^?^yi( r i' r 2? r 3^ "t = r i — r 2), i-G., the density probability of finding a particle of 
species % at r 3 in the presence of the dumbbell. Mathematically the conditional three particle 
distribution function, ^^(r 3 ;t), is related to the homogeneous three particle distribution 
function fi^ri, r 2 , r 3 ) by 
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(3) I \ 

W«t) = -^— P) 

[31 (2) 

The projection of <7kL( r i> r 2 , r 3 ) gives directly gf^j((t). This projection can be provided by 
the Born- Green- Y von theorem (BGY) that is based on a balance of the mean effective force 
F^[^(r 3 ;t)]. 

B. The Born- Green- Yvon equation or a force law 

The Born-Green- Yvon (BGY) equation is one of the so called hierarchy equations and it 

nn 

is an exact theorem relating the n and (n+ 1) particle distribution functions [£9j,|ZQj]. Here, 
we derive the BGY equation as a sum of all the forces exerted on one of the two dumbbell's 
particles (let us say particle of species 7 at r 2 ). This mean force has two contributions: 
(i) the direct force f/3 7 (t) exerted by the particle of species (3 at r 2 and (ii) the force f 7 (t) 
exerted by all the other particles. Thereby, the total mean force F^ 7 (t) reads 

F^ 7 (t) = f /37 (t) + f 7 d (t), (8) 

Assuming that the dumbbell and fluid species are spherical particles interacting through 
central force potentials, the component of fg 7 along t is given by 

tM r = IM) = ^ — , (9) 

where u^{t) is the potential of direct interaction between the two dumbbell parti- 
cles. The elementary force df^ produced by a fluid element at r3 is given by df~ = 
I2i=i f-yi( r 23) Pi ( r 3) dv 3, where f 7 i(r 23 ) is the force between a particle of species i [of local 
density Pi(r 3 ) = Pigfl fi (r 3 ; r)] at r 3 and the dumbbell's test particle of species 7. The 
component of dfz along the direction of t is given by 

< = E ^ ( r 23) Pi (rs) dv 3 = - £ t-r 23 ^J (r23) p t (r 3 ) dv 3 , (10) 
i=i i=i z ° 

with t and r 23 being unit vectors along the t and r 23 directions, respectively, w 7 i(r 23 ) is the 

potential of interaction between an z-species particle with the 7-species particle. Substituting 



6 



Eqs. © and JUJ into Eq.@, F^ 7 is given by [33j,|37| 



/ \ dwp 7 (r) dup 7 (r) A f du 7 j(r 23 ) [3] 



5>/ ^^cos^ 3 ;(r 3 ;r)^3, (11) 



dr dr f=x ■> dr 23 

where t-r 23 = cosQ and w^{t) is the potential of mean force between the two dumbbell's 
particles. According to Eq. ((H), 



and thus 



w P~f( T ) = ~ k BT In [gpj(r)\ , (12) 



dinger) du Pl {r) * r du yi {r 23 ) [3] 
W^-— -— ^— -X>y cos fl^ for)**, (13) 



which is the Born-Green- Y von (BGY) equation. The degree of accuracy of <?/3 7 (t) depends 
on the method used to compute g l oL (r 3 ;r). If Ogi^ (r 3 ;r) is computed through the TPE 
of integral equation theories, it was shown that new diagrams are included in the cluster 
expansion of g^lir) 3?1 (see Fig. |3J). By examination of the transformation of the Mayer 
diagrams through the formalism outlined above, the denomination of three point extension 
becomes clear. A more detailed description of TPE can be found in ref. [37]. 

C. Application to the RPM electrolyte 

In the RPM electrolyte the fluid is considered as made up of hard spheres of diameter 
a with a central charge g« = ^e, where is the valence of species i and e is the protonic 
charge. The electroneutrality condition for the n-component electrolyte is 

n 

£ w = °- ( 14 ) 
1=1 

Assuming that the dumbbell particle (a species) is made up of two particles of the same 
species from that in the fluid (see Fig. |2J), the TPE-HNC/MSA equations become 

(r u )dr 4 L (15) 

where 

u ai (r 3 ) = u ai (ri 3 ,r 23 ) 



m=l 



7 



2 2 

1 '■ — it r 13 and r 23 > a 



OC 



£r 13 



(16) 



if n 3 or r 23 < a 



with zp and z 7 standing for valence number of particles (3 and 7, respectively. For spherical 
ions the direct correlation function depends only of the ions distance r 34 = |r 34 | . Within the 
mean spherical approximation, its analytical expression is 



Crm(r 34 ) = c hs (r 34 ) + z m ZiC sr (r u ) - (3 



(17) 



where c hs (r 34 ) is the direct correlation function for a hard spheres fluid in the PY approxima- 
tion and c sr (r 34 ) is a short ranged function. Because of the s ymm etry around the dumbbell 
axis, it is convenient to use prolate coordinates (77, £, <fi) 0,lzi] defined as follows 



x 

y 

z 



2 

1* 



(7? 2 -l)(l-aOOS0, 

(r/ 2 -l)(l-asin0, 



(18) 



and where the volume element is given by 



dv = — (rj 2 — £ 2 ) dcf)d£dr). 
8 ^ ' 

The relative distance r 34 is then given by 



(19) 



rl = j { (nl - 1) (1 - it) + (vl - 1) (1 - il) + - 

- 2y/(rj£ - 1) (1 - 0) (Vl " 1) (1 - il) cos0 4 } . 



(20) 



In prolate coordinates the potential of electrostatic interaction, iQ, between the dumbbell 
and one fluid ion of species % can be conveniently rewritten as 

and Eq. ([To]) as 
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( t \ [3] / c \ J 2 P e2 ( Z P Z ~1 

g ai (773, 63) = 9/3ji (vs, 63; r ) = ex P \ — — — —r + ~ ~ 



t£ \m -63 m + &J 

rl rco 

/ / Pas (774, £4) K (773, £ 3 , 774, £3) dr7 4 d£ 4 

•7-1 " / '?0(?4) 
/■l poo 

/ / Pad (7/4, 64) L (773, 6, 7/4, &) d77 4 d^ 4 

rl roo ^ 

~ Zi / p Q d (774,^4) A (773, 63,774, C4)d774d^ 4 - J (773,6) h (22) 



+ 
+ I. 



with 

+ for 6o<6<1 
1 for < 6 < 60 

and 60 = 1 — b = 2a/r and 77o(— 6) = 770(6)- The expressions for K, L, A, J, p as and p a d 
are 

K(77 3 , 6s, 774, 6) = T (t7 4 2 - 6 4 2 ) / c hs (r 34 )# 4 , 
8 jo 

7"^ /'(/'max 

1^(773,6,774,64) = -^-(774-64) y o c sr (r 34 )# 4 , 



A 773,63,774,64) = 5 — (774-64) / — . 

Pas (774, 64) = P/3 7 s (774, 64) = PrnKm (^4, 64) , 

m=l 

n 

Pad (774, 64) = P/3 7 d (774, 64) = Z mPmh am (t]A, , 



4 



m=l 



J (773, 63)=/ / K(77 4 ,64,773,63)d77 4 d64+ / / K (774, 64, 773, 63) d77 4 d64, 

respectively, with 

6min(r) = < 



if r < a 
60 if t > a- 

By introducing the elliptic function of second kind F(tt/2, k), one can rewrite A as 



A(r, t 77 n - T {VI - ® F (?r/2 ' k) (23) 

A (773, 43, 774,447 - {M) 



^'34 

where 



fc2 = rV(77 3 2 -i)(i- 63 2 ) (vl - 1) (i - 61) 

n/ ^max\9. ' \ / 



2(r 3 m 4 < 



9 



and 

2 + - V±U) 2 • (25) 

Eq. EH is in fact a set of two coupled, three dimensional, non-linear integral equations. To 
solve these equations, we have developed a sophisticated, but efficient, finite element method 
for its solution (see appendix for details on our numerical method). 

Using Eq. (fTTjl . the mean force between the two dumbbells particles reads 

Ffh{r) = fa (r) + $, (r) , (26) 

where 



( maxy 
I' 34 I 



T 

T 



V^i - 1) a - a) + v^-iMi-a 



#r W = E Pi /' , A <ft ^ &) > & r] [-2£f - 3fo£ 3 2 + (2 - b") £ 3 + 6] dfa (27) 



and 



M T ) = li + — : — / . / .^Mftiw - — — ^3^ 3 - (28) 



Thus the pair distribution function of the electrolyte solution is given by 

9fh {r) = exp {-/3 £ F^(r)dr} . (29) 

The solution of Eq. (j22)) and calculation of through Eqs. ()26)). (|77jl and (J28J) were 
numerically solved. 



III. MOLECULAR DYNAMICS 



The electrolyte is confined in a cubic box of length L. The bulk salt concentration p is 

TV 

then given by — , where N is the number of positive (or negative) ions. The dumbbell is 
made up of two fixed ions (with a center-center separation r) disposed symmetrically along 



the axis passing by the two centers of opposite faces. A similar system setup was also used 
elsewhere to study two fixed macroions [3, [z2, [z^J. We use MD simulations to compute 
the motion of the mobile fluid ions coupled to a heat bath acting through a weak stochastic 
force Wj(i) with a zero mean value. The equation of motion of any mobile ion i reads 
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m^ = -V t U-mT^ + W l (t), 
at 2 at 



(30) 



where m is the ion mass, T is the friction coefficient and — ViU is the potential force 
having two contributions: (i) the Coulomb interaction and (ii) the excluded volume in- 
teraction. Friction and stochastic force are linked by the dissipation-fluctuation theorem 
(W f -(t) • W,(t')) = QmTk B T6 l3 6 (t - t'). 

Excluded volume interactions are modeled by a pure repulsive Leonard- Jones (LJ) po- 
tential defined by 



U LJ {r) = 



4exj 
0. 
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(31) 



+ e LJ , for r < 2 1/6 a 

for r > 2 1 / 6 o, 
where a is the ion diameter. 

The electrostatic interaction between any pair ij, where i and j denote either a dumbbell 
ion and/or a mobile fluid ion, reads 



U d {r) = ±k B T£ B -, 



where +(-) applies to ions likely (oppositely) charged, £ B 



ek B T 



(32) 



is the Bjerrum length 
- Zj = z). To link our 



describing the electrostatic strength and z is the salt valence (zj = Zj = 
system parameters to experimental units we choose the LJ energy parameter e^j = k B T 
(where T = 298K) and a = 4.25 A. This leads then to the water Bjerrum length £ B = 1.68a = 
7.14 A. A macroscopic system was mimicked by imposing periodic boundary conditions. 
The long range Coulomb interaction was treated by using an optimized and efficient Ewald 
summation variant, namely the particle-particle-particle-mesh (P3M) method [3]. 

In order to limit the size effects, we choose L sufficiently large, typically 10 times (or 
more) the Debye-Hiickel screening length. The number of ions in the box is 500 for all cases 
(concentrated and dilute solutions). It is important to mention that the computation of 
gpL(r, 9',t) is statistically extremely demanding and especially for small 9 angles, since the 
quantity of information varies like sin(#). In this notation, the distance r = \ts\ and the 
angle 9 = Z(r!,r 3 ) are always relative to the center of the dumbbell (see Fig. |2J). The fact 

[31 

that the observable gpL(r, 9; r) concerns only an "elementary solid angle" , it strongly reduces 
the available information compared to that available for the pair correlation function, since 
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in that latter case a full solid angle 4ir and many ion pairs are accessible. To overcome 
this difficulty, we considered a sufficiently large angle range A9 (typically between 5 — 15° 
depending on the concentration p), so that the gathered informations contains as less noise 
as possible. On the other hand, A9 must not to be too large otherwise the resolution gets 
too small. For each system under consideration, a compromise between these two effects 
that had to be found. 

Finally for the computation of the effective mean force between two ions, we considered 
the same system but where no fixed dumbbell is present. Thereby, we could compute the 
potential of mean force, knowing the g(r), and then get by derivation the effective force. 

IV. RESULTS 

We have done calculations for the 1:1 and 2:2 electrolytes using the TPE-HNC/MSA 
integral equation. For the 1:1 electrolyte the agreement between TPE-HNC/MSA and MD 
results is qualitative and quantitatively very good. However, in order to keep low the number 
of plots we just present a detailed analysis on the results of the 2:2 electrolyte. The choice 
of divalent ions is motivated by the fact that it represent a strong test for liquid theories. 
Thereby, we considered two typical concentrations: (i) the concentrated case with p = 1M 
and (ii) the dilute case with p = 0.005M. As a main result, the effective mean force obtained 
by TPE-HNC/MSA, HNC/MSA, and MD simulation is presented for each concentration 
regime. In order to further quantify the robustness of the TPE-HNC/MSA theory, we 
investigated in detail the conditional three-particle distribution function, g^^r, 9;t), by 
comparing TPE-HNC/MSA with MD. 

For the discussion, it is convenient to adopt the following notations: g + \ (r, 9; r) stands 

for the distribution function of negative ions when the dumbbell is made up of two positive 
[si 

ions, g + (r,9;r) for that of negative ions when the dumbbell is made up of a negative 

and a positive ions, and so on. By symmetry the three particle distribution function sat- 
isfies <7++_(r, 9; r) = g^_ + (r, 9; r) and also (r,9;r) = #+L + (r, 7r — 9;t). Thereby, we 

systematically compared theory and simulation for gf^r, 9; r), but show results only for 
9p},i(r, 9;t = a) (i.e., when the two dumbbell ions are in contact), for two given values of 
9 (7r/4 and 7r/2). In addition, within the TPE-HNC/MSA theory, we also provide the full 
^-dependence of gf\i(r, 9; r) for different r. 
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A. Concentrated case 



In this section, we deal with the concentrated electrolyte solution (p = 1 M). The elec- 
trostatic screening at such high ionic density and valence {z — 2 ) is very strong. The study 
of such a system is important to test TPE-HNC/MSA theory, since already inhomogeneous 
and homogeneous HNC/MSA theories are in excellent agreement with molecular simulations 



under such conditions 



76|. 



1. Three particle correlation function 

a. Symmetric dumbbell We first consider symmetric dumbbells made of like charged 
positive divalent ions. The profiles of g++_(r, 9 = tt/2; t = a) and g++ + (r, 9 = n/2; r = a) 
are depicted in Fig. 0J Concerning the negatively charged fluid ions (i.e., "dumbbell counter- 
ions"), we have quantitative agreement between theory and simulation even near the distance 
of closest approach. The slight difference at short distance (r ~ a) is due to the fact that 
for the short-ranged excluded volume interaction, MD simulation is built with a soft-core 
LJ potential whereas the actual theory uses a true hard-core potential. For the positively 
charged fluid species ("dumbbell co-ions"), we also have an excellent qualitative agreement. 
The TPE-HNC/MSA maximum of the co-ion distribution function is within the statistical 
error, however slightly higher than MD data. The location of the maximum is nearly the 
same as that found with simulation. Hence TPE-HNC/MSA has an excellent agreement 
with MD, within the numerical error. 

For 9 = 7r/4 (see Fig. [SJ, one still has the same quantitative agreement between 
MD and TPE-HNC/MSA for the dumbbell counter-ions. It is observed that the value of 

g + \ (r, 7r/4; a) at closest approach (r = 1.29a) is not as high as at 9 = it/2 (see Fig. EJ) for 

the corresponding plot. The physical reason of this feature is straightforward. The closest 
approach to the center of the dumbbell is larger at 9 = tt/4 than at 9 = tt/2, therefore, since 
all particles have the same size, the resulting attractive electrostatic interaction between the 
dumbbell and the counter-ion is higher at 9 = it/ 2. For the dumbbell-co-ions distribution 
<7+!|. + (r, 7r/4; a) we have quantitative agreement between TPE-HNC/MSA and MD. 

The three dimensional (3D) plots of the three particle (counter-ion-dumbbell) distribution 
function g [ + ] + _(r,9;r) obtained by TPE-HNC/MSA are sketched in Fig. El At r = a 



13 



(dumbbell ions at contact), Fig. EJ^a) shows a strong variation near to the surface of closest 

approach. As expected, the maximum is obtained at 9 = 7r/2,37r/2 (g + \ « 50), whereas 

the minimum is at 9 = 0, it (g+\ « 8). Moreover, we have oscillations in the distribution 

function, as a function of r, for any 9 , which confirms our previous observations at 9 = 
7r/2,7r/4 (see Figs. EJandEJ). We have carefully checked that this feature holds at any r. 

roi 

The 3D plot of the co-ion-dumbbell distribution g +++ (r, 9; a) is not reported here. 

At a larger dumbbell separation, r = 2a [see Fig. E](b)], g++ + (r, 9; 2a) is still highly peak 
at 9 = tt/2 and has its maximal value at the middle point of the dumbbell. At sufficiently 
large separation, we have an isotropic counter-ion distribution around each dumbbell particle 
(not shown here). We point out that although the probability of finding two like-charged 
ions in contact ( r = a) is very low, the probability of having more than two counter-ions in 
contact (at 9 « n/2) with the two like-charged ions dumbbell, is very high. This implies an 
overcompensation of the dumbbell's charge, which is verified by the observed oscillations in 
the counter-ions profile of Fig® since oscillations imply an electrical field inversion, which 
implies charge reversal. It should be stressed that to calculate thermodynamics functions 

roi 

such as the internal energy or pressure, g ++i (r, 9; r) for every r must be known, even if the 
probability of finding two like-charged ions at contact is very low. 

b. Antisymmetric dumbbell We now consider antisymmetric dumbbells made of two 

opposite divalent ions. In this case, by symmetry arguments we expect that (r, tt/2; a) = 

g+L + (r, 7r/2; a). The profiles of g+. (r,n/2;a) and #+L + (r, 7r/2; a) are plotted in Fig. 

Since at 9 = n/2 the electric field component (produced by the dumbbell) perpendicular 
to the dumbbell axis is zero, the electrostatic correlations are only generated by the fluid 
ions. Consequently, we expect a quasi-neutral fluid behavior. This is precisely what Fig. [7| 

shows for theory and simulation, where g + (r, tt/2; a) and g + _ + (r, 7r/2; a) curves collapse 

in a single curve. The adsorption at contact, is a hard sphere entropic effect due to the salt 
high concentration. This adsorption does not occur at low salt concentration. 

Results for 9 = tt/A are shown in Fig. El We have again a very satisfactory agreement 
between theory and simulation. 

The 3D plot of #JL(r, 9;t = a) = g [ +l+(r, n-9;r = a) obtained by TPE-HNC/MSA 
is sketched in Fig. El The maximum and minimum are located at 9 = and 9 = n 
at dumbbell contact, which implies a high probability of a line quadruplet configuration. 
However, if we look at FigJ3 it implies that, although with a lower probability, positive or 
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negative ions are adsorbed around the center of the antisymmetric dumbbell, at 9 = 7r/2. 
Again, we observe oscillations at any 6 angle, and we checked that it is the case for any r. 
Hence, charge reversal is also present, implying that more than two ions are adsorbed to the 
dumbbell. Thus, probably compact clusters more than line clusters are formed at this high 
concentration. We will come back to this point later. 



2. Effective force 

The effective mean force between two like charges [F ++ (r)] and that between two opposite 

charges [F + _(r)] as a function of their mutual separation r are depicted in Fig. [TU1 in reduced 
kjfT 

units of — — . As expected, theories (TPE-HNC/MSA and HNC/MSA) and simulation are 

"b 

in very good agreement for both forces F ++ (r) and _F + _(r). 

An interesting feature is the kink in F ++ (r) occurring at r = 2a, that is somewhat less 
marked, however present, on the simulation plot (due to the softness of the ions and also the 
lower radial resolution there). This jump in the first derivative, F+ + (r), at r = 2a is not an 
artifact of the theory (or the simulation) but a true physical feature. This effect is due to 
excluded volume correlations and, in much lesser degree, to electrostatic correlations. It is 
clear that, at r = 2a, the configuration consisting of a counter-ion lying exactly between two 

co-ions (i.e., H h) is energetically very favorable (see Fig. Eh). This implies the formation 

of ion complexes, in qualitative agreement with Caillol and Weiss [l^J and Yan and de Pablo 



2l|. When r > 2a (more precisely r — > 2a + ), the presence of an in-between counter-ion 
leads to a relatively strong resistance, on the level of the depletion force, upon approaching 
the two co- ions. On the other hand, when r < 2a (more precisely r — > 2a~ ), the absence 
of an in-between counterion leads to an easier approach (on the level of the depletion force) 
of the two co- ions. These mechanisms, explain (i) the discontinuity of F' ++ (r) at 2a and (ii) 
the fact that \F+, (r — > 2a~)| < \FL,(r — > 2a + )|. This effect should also be observed in 
neutral hard spheres systems at sufficiently high density, and in the interaction between two 
macroions. 

As far as the force F + _(r) is concerned, this kind of discontinuity in the derivative is 
absent or nearly undetectable. This is due to the fact that, at r = 2a, the probability of 

finding the configuration consisting of an ion between two oppositely ions (i.e., H ) is 

considerably smaller compared to that obtained with the configuration H h -F+_(r) < 
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and F ++ (r) > are of the same order of magnitude and indicate, of course, that the (H — ) 
configuration is of high probability, whereas the (++) is of low probability. 

By definition In gpL(r, 9; r) = — w^i(r, 9, r)/k B T. Hence, PigfL(r,9;T) gives the proba- 
bility of finding an ion of species i, at a certain position (r,9), from a dumbbell made of 
two ions of species j3 and 7, located at a distance r, from each other. wp y i(r, 9,t) is the 
potential of mean force between the ion % and the dumbbell. The mean internal energy of an 
ideal gas per particle is 3k B T '/2. Hence, —Wo = — w ++ -(r, 9, t)/[7cbT] > 3/2, for a plus-plus 
pair, i.e. g+\__(r, 9; r) > 4.48, implies that a counter-ion next to a like-charged dumbbell has 
an adsorption energy larger that its thermal energy, and thus it is tightly attached to the 
dumbbell. For a 1M electrolyte, -W > 3/2 implies p_g [ + ] + _(r, 9; r) > 4.48 M. In FigEk, 
the peak is for p_g l + ] + _(r ,9 = vr/2;r = a) = 50M > 4.48M, i.e., -w ++ _(r,9 = tt/2,t = 
aj/ksT = 3.9 > 1.5. On the other hand, at 9 = ir, p-g ++ _(r , 9 = n; r = a) =5M, i.e., 

— u? +H (r, 9 = tt,t = a) /k B T— 1.6~1.5. Therefore, a negative ion will be strongly attached 

to the positive ions pair (at 9 — tt/2). A simple calculation shows that the unscreened 
attractive electrostatic energy of a second negative ion to the (+H — ) ion complex decreases 
to around 50% of the attractive energy of the positive ion pair to the first negative ion. 
Hence, a second adsorbed ion, at 9 = ir/2, seems likely. Thus, Fig. suggest a quadruplet 
structure, where the two counter-ions are at 9 = tt/2. Clearly, more than two counter-ions 
are adsorbed, since the dumbbell charge is overcompensated, i.e., there are concentration 
profile oscillations. The adsorption of these additional counter-ions is due to the short range 
correlations, i.e., ions next to the dumbbell feel a net force toward it due to the uneven 
collisions from bulk ions, and is an entropic effect, beyond the ideal gas entropy. This effect 
is larger, the larger the electrolyte concentration, and it will not be present in a point ion 
electrolyte. Because the attractive potential of mean force is very high, this compact ion 
complex structure is very stable, although very unlikely, because F ++ (r) > (see Fig. ITTH). 

However, for 2a < t < 3a, a (H h) configuration is very likely. Hence, this indicates that 

there are several mechanisms for the formation of ion complexes. 

For the plus-minus pair, in Fig [HI the peak is for p+g + (^0?^ = 0;r = a) = 6.8M > 

4.48M, i.e., -w + __(r,9 = 0, r = a)/k B T = 1.9 > 1.5. and from Fig. □ -w + __(r,9 = 
tt/2, t = a)/ k B T = 0.18 <§; 1.5. Hence, for an unlike charged dumbbell, we expect an aligned 
stable quadruplet configuration (because of symmetry), due to energy arguments. However, 
due to entropic effects more counter-ions are adsorbed into the dumbbell, producing charge 
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reversal, as can be seen from the oscillations of FigEl These additional ions are delocalized 
around the dumbbell, hence, generating compact ion complexes because F + _(r) < (see 
Fig, imp, this configuration is very likely. 



B. Dilute case 



In this section, we study a dilute, divalent electrolyte (p = 0.005M and z = 2). To the 
best of our knowledge all of the known liquid theories fail to describe the RPM behavior 
under these conditions [?], Hence, the study of low concentrated solutions of multivalent 
ions represents a strong test case for a liquid theory. In addition, for the RPM electrolyte 
we are on the low concentration side of the phase diagram. 



1. Three particle correlation function 

a. Symmetric dumbbell Figures UTT a) and [TTT b) show a comparison between TPE- 

HNC/MSA and MD results for g+\ (r, ir/2;a) and g^ ++ (r, 7r/2; a), respectively. One can 

see that the electrical double layer is wider than in the concentrated case i.e., the correlations 
are long ranged. For g+ + -(r, vr/2; a) [see Fig. ITTTa)]. TPE-HNC/MSA and MD results show 
quantitative agreement, even near contact. For the g++ + (r, 9 = 7i/2;a) [see Fig. ITTT b)]. a 
qualitative agreement between TPE-HNC/MSA and MD results is found. 

Pol 

At 9 = it/A (see Fig. H2J) it is found that the contact value of g ++ _{r, 7r/4; a) (about 500) 
is much smaller, of two orders of magnitude, than that at 9 = tt/2. This can be explained 
in terms of the electric field produced by the dumbbell at 9 = it/ 2 which is considerably 
stronger than at 9 = 7r/4. 

The 3D plot of ln(g+\__(r, 9; r)) can be found in Fig. [T31 For r = a [see Fig|T3] (a)], 
it is observed a strong variation of the distribution function close to the dumbbell (at the 

surface of closest approach). As expected, the maximum of g+\ (r, 9; a) is at 9 = tt/2. On 

the other hand, at r = 5a (see Fig. H~3T b)). the angular variation of g+\ (r, 9; 5a) (near 

contact) around one ion of the dumbbell is not as peaked as in g + \ (r, 9; a). However, the 

dumbbell ions are still correlated, i.e., their electrical double layers are strongly overlapped 
although r = 5a. We had to go up to r = 60a (not shown) to cancel the overlapping of 
the electrical double layers of the dumbbell ions. At low salt concentration there is a longer 
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range penetration of the ions electrical field into the fluid, and hence charge correlations 
are of longer range. For low salt concentration the role of higher order diagrams is more 
important. It is observed that for this low concentration case, there are no oscillations in 
the counterion concentration profiles. Hence, no charge reversal is present and, thus, the 
formation of a simple more complex ionic configurations, beyond a quadruplet formation, is 
not supported by our results. Because of the very large value of Wo ~ 8.6, the adsorption 
of two counter-ions to the like-charged dumbbell (at 9 = tt/2) is much larger than for the 
equivalent situation for the concentrated case, where W ~ 3.9. Hence, for the dilute case 
the quadruplet is more stable, but even less probable due to the lower concentration. 

b. Antisymmetric dumbbell We now consider the three particle distribution function 
where the dumbbell is made up of two opposite ions, and for the same fluid parameters 
as in Figs. fTTI and fT2l Only the case of 9 = tt/A is shown (see Fig. 114*)) . given that for 
9 = 7r/2 the electrical field is zero and since the electrolyte concentration is very low we 

have (r, 9 = tt/2) = g+L + (r, 9 = vr/2) « 1. This is in contrast with the 1M electrolyte 

result of Fig. For 9 = n/2 the same good agreement is found between TPE-HNC/MSA 
and MD as that in Fig. H3 

The 3D plot of ln<7+! (r, 9\r = a) = In g+L + (r, n — 9; r = a) which is the potential of 

mean force is sketched in Fig. This function is quasi center-symmetric with respect to 
the dumbbell center. This feature is due to (i) the symmetry of the electrostatic correlations 
and (ii) the fact that the contribution of the excluded volume correlations (at such low 
density) is negligible compared to that in the concentrated case. As expected the function 
is strongly peaked at 9 = 0. For the concentrated case (not shown) the asymmetry is 

higher. The important result shown in this figure is the formation of a stronger (H 1 — ) 

line quadruplet, than for the concentrated case, since here — (r = 3/2a,9 = 0, r = 

a)/ksT = — w + _ + (r = 3/2a,9 = ir, r = a)//c^T w 4.2 > 1.5, which is much higher than 

that for the corresponding concentrated case (— u> H (r = 3/2a,9 = 0,r = a)/kBT w 

1.6). Hence, this line quadruplet structure would be more stable. This result suggests that 



quadruplets, if present, would be in a linear configuration more t 
structures, in disagreement with Fig. 4 of Yan and de Pablo 



ran in compact quadruplets 
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2. Effective force 



The effective mean force between two like charges [F++(r)] and that between two opposite 
charges [_F + _(r)] as a function of their mutual separation r can be found in Fig. in 
reduced units of Concerning F + _(r), theories (TPE-HNC/MSA and HNC/MSA) and 

simulation are in quantitative agreement. 

As pointed out above, for F ++ (r) in the concentrated case, the derivative F+ + (r) is again 
discontinuous at r = 2a. The same mechanism proposed for the concentrated case (see 
Sec. IIV A 2|) applies here. This important feature is not captured by HNC/MSA, proving 
the qualitative improvement by using the TPE method. This better description steams from 
proper inclusion of long ranged correlations. Finally, we have a good quantitative agreement 
between TPE-HNC/MSA and MD. In comparison of Fig. with that for the concentrated 
case, Fig. two important differences are observed: The intensity and the range of the 
force is larger for the dilute case, implying that the electrical field is less screened. Also 
for low concentration F ++ > 0, i. e., it is always repulsive, whereas in the concentrated 
regime, for some interval of r, F ++ is negative, implying an attraction and hence different 
ion-complexes mechanisms. In addition, one can expect that for a certain combination of 
temperature, solvent dielectric constant, and salt valence and low concentration, one can 
find a phase transition, in which associated ions and free ions coexist: single ions, ion pairs 
and quadruplets. Hence, from Figs. 01 and for the concentrated case, and Figs. E3 and 
[HE for the dilute case, we conclude that linear ion complexes are likely to be formed. At 

low concentration, dumbbells (H — ) and line quadruplets (H 1 — ) are very likely to be 

formed, whereas at high concentration larger complexes than quadru plet s are formed. This 
is in qualitative agreement with the predictions of Caillol and Weiss jlTf. 



V. CONCLUSIONS 



We have investigated the structure of 1:1 and 2:2 RPM electrolytes by means of integral 
equations and MD simulations. Using the three point extension to the HNC/MSA theory, 
the conditional three particle distribution function, g^(r, 9; r), was computed and compared 
with that obtained by MD. Although it is not shown, for the 1:1 electrolyte the quantitative 
agreement between TPE-HNC/MSA and MD is excellent. For the 2:2 electrolyte, we explic- 
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itly report here results for two typical concentrations: (i) the concentrated case (p = 1M) 
and (ii) the dilute case (p = 0.005M). 

As far as the concentrated case concerns, it was found that g' 3 l(r, #) always presents 
oscillations. The detailed comparison between TPE-HNC/MSA and MD, carried at fixed 
separation r = a (between the two constitutive ions of the dumbbell), shows an excellent 
qualitative and/or quantitative agreement. This is true for all values of r (not shown). 

On the level of the effective mean force between two ions, both, TPE-HNC/MSA and 
HNC/MSA are in very good agreement with MD. This is consistent with previous compar- 
isons between HNC/MSA and Monte Carlo results . Hence we can conclude that the 
TPE-HNC/MSA method is also suitable to describe concentrated electrolyte solutions. It is 
important to point out a particular behavior in the effective force between like-charged ions 
[F ++ (r)] observed at r = 2a, where an abrupt change in its slope appears due to excluded 
volume correlations. This behavior can not be directly seen in the pair distribution function 
(for this value of p). 

In the dilute regime, the analysis of the three particle distribution function and the effec- 
tive force shows the long range nature of the correlations. For the three particle distribution 
functions, we had to go up to a distance separation of r « 60a, in order to uncorrected 
the two constitutive dumbbell ions. Again a good agreement for g^(r,9) is found between 
TPE-HNC/MSA and MD, proving the robustness of the TPE formalism. The study of the 
effective force reveals a quantitative agreement for the force between two oppositely charged 
ions, F+_(r), between TPE-HNC/MSA and MD, although HNC/MSA is also very good. 
For the force F ++ (r) we again remark the occurrence of an abrupt change in its slope at 
r = 2a, which is not predicted by HNC/MSA. On the other hand, TPE-HNC/MSA and MD 
are in quantitative agreement, showing the ability of TPE to take fairly well into account 
long range correlations. It is precisely in this region of the ionic fluid phase diagram, i.e., 
low concentration and high Coulombic coupling, where all the other theories fail. 

The TPE approach is a general formalism that improves existing liquid theories, by 
including higher order diagrams in a systematic, consistent way Here we have shown it 
to be successful for ionic fluids, in all the regions of the RPM phase diagram, in particular 
in the region of low salt concentration and high coulombic coupling. 

In the high concentration regime, ion pairs tend to form aligned quadruplets, i.e., ( — h 
— h) structures are energetically favored. However, because of short range correlations, 
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other delocalized ions are adsorbed to produce charge reversal of the unsymmetrical ion 
dumbbell and hence the formation of larger complexes than quadruplets is favored. In the 

low concentration regime, the ( — I h) aligned quadruplet structure is even more stable 

than for the high concentration case. Hence, dumbbells and aligned quadruplets are likely 
to be formed. No adsorption of additional ions is present, since there are no oscillations in 
the concentration profile and, hence, there is no charge reversal of the dumbbell or higher 
multiploles. In the high concentration regime charge reversal is present, whereas at low 
concentration there is no charge reversal. Our results clearly indicate the formation o 



pairs and complexes, in agreement with previous theoretical predictions 
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581 ] and simulation results 
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e rormation oi ion 



24j . In our theory we do not impose ion pairs, and 
could be useful to explore RPM phase transitions, critical behavior and could provide a 
means to understand the molecular mechanisms behind fluids phase transitions. 
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APPENDIX A: NUMERICAL METHOD 
1. Finite Element Method 

To solve the TPE-HNC/MSA equation, Eq. (|22p. it is necessary to use a numerical 



method, since an analytical solution is not available. The finite element method 
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EM) has 
and it 



been used in the past to solve HNC/MSA equation in several geometries {35 
has been proved to be efficient. The general form of TPE-HNC/MSA integral equation can 
be written as 



g ai (rj, = exp M f (77, + / [°° E PmKm (rf, £') F ( V , f , rj', £') drfd? , (Al) 

(A2) 
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where h am (r},£) and Mi(rj,£) are functions defined on a bidimensional domain (77, £) G 
[—1,1] x [l,oo). Since h am (r),£) 7^ only in a region close to the dumbbell, we solve 
Eq. (|A1J) just in a finite domain. In the FEM [3| , the domain is divided into N elements. 
Every element in a domain Ax is divided into L sub-elements. In prolate coordinates, the 
dumbbell geometry of Fig. El is mapped into the geometry shown in Fig. EE where one of 
the N triangular elements is shown. 

In order to solve Eq. (|A1|L the function h am (77, £) is expanded as a linear combination 
of a L base elements |0f (77,0 ,7 = 1, ...,L \, where 1 < K < N. These base functions 
are defined in such a way that 4>f (77, £) = if (77, £) ^ Aff. Furthermore the base functions 
are chosen so that for a set of L points (rjj,^j) (which are called nodes, see Fig. ITS)) , they 
satisfy 

0fO7sO = ( %> with i,j = 1,...,L -, (A3) 
with 5jj being the Kronecker delta function. Hence, 

JV L 

h am (r/, = EE "*i<f>f (V, ,m = 1, 2 (A4) 
x=i z=i 

where |o;^,/ = l,...,Lo| are the L$ coefficients of the h am in the i^-th finite element. 
Thereby, the coefficient is the value of the function at the i -th node, i.e., 

^mi = ham(Vh€i) ■ ( A5 ) 

It is useful to renumber (pi (77, £) and ui m i, so that Eq. (|A4|) can be rewritten as follows 

Mo 

Km (T/, = E fa > 111 = !> 2 ( A6 ) 

2=1 

where 

w mi = ^am (^i ; £i) = 3am fa, £i) — 1) ( A 7) 

with M = N x L , and the superscript X has been omitted. By substituting Eq. ()A6|) into 
Eq. (IA"T| . we get 

{Mo 2 ] 
Mi (77, + E E Pm^mlCl(77, , (A8) 
1=1 m=l J 

where 
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Q (17, = / J <h (jf, O ^ (»/, »?, «W (A9) 
Evaluating Eq. ()A8|) at the fc-th node an d using Eq. (jA7|) . we get 

= exp < M ik + Pm^miQk \ ~ 1, (A10) 

with Cjfc = Ci(r)k,£k) an d Mjfc =M* (77*., . Thus, we have a system of 2M non-linear 
algebraic equations which can be solved by any of the standard methods, for example the 
Newton's method, which is our method of choice. 



2. Choice of the base functions <pi (77, £) 

To construct the base functions, it is necessary to use a coordinate system defined in the 
element's domain. In our method, we have used the area coordinates (Li) defined as follows 

a,- + bill + c,£ , . s 

Li = ^ with 2 = 1,2, 3, (All) 

where 2A is the area of the triangular element and 



"i = '///■• - '//. s. r 

&i = &-&, (A12) 
Q = Vk - Vji 

with cyclic rotation of indexes, where j, k = 1,2,3 but i 7^ j 7^ k. The set of points 
{(Ci,Vi) , * = 1,2,3} are the coordinates of the triangle corners. The relation between the 
coordinates (£, 77) and the triangular coordinates {Li, i — 1,2, 3} is given by 



T] = Lit]! + L 2 7] 2 + L3773, 

£ = L^ + L^ + Lsts, (A13) 
1 = L1 + L2 + L3. 

The number of nodes is equal to the number of base elements. A quadratic base was used 
to solve Eq. (jAl|) and therefore L = 6. For the corner nodes we have 
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4>i = (2Lx — 1) Li, etc., 



(A14) 



and for the middle- side nodes 



04 = 4LiL 2 , etc. (A15) 
In this coordinate system, Eq. (jA9|) becomes 



/o Jo 

xF(L l ,L 2 ,r ] ,£)dL l dL 2 . (A16) 



[1] J. C. Rasaiah and H. L. Fiedman, J. Chem. Phys. 50, 3965 (1969). 

[2] J. C. Rasaiah, Chem. Phys. Lett. 7, 260 (1970). 

[3] J. C. Rasaiah, J. Chem. Phys. 52, 704 (1970). 

[4] J. C. Rasaiah, J. Chem. Phys. 56, 3071 (1972). 

[5] E. Waisman and J. L. Lebowitz, J. Chem. Phys. 56, 3086 (1972). 

[6] E. Waisman and J. L. Lebowitz, J. Chem. Phys. 56, 3093 (1972). 

[7] D. Henderson, M. Lozada-Cassou, and L. Blum, J. Chem. Phys. 79, 3055 (1983). 

[8] T. L. Croxton and D. A. McQuarrie, J. Phys. Chem. 83, 1840 (1979). 

[9] D. N. C. J. C. Rasaiah and J. P. Valleu, J. Chem. Phys. 56, 248 (1972). 
[10] D. N. Card and J. P. Valleau, J. Chem. Phys. 52, 6232 (1970). 
[11] S. A. Rogde and B. Hafskjold, Mol. Phys. 48, 1241 (1983). 
[12] D. A. McQuarrie, J. Phys. Chem. 66, 1508 (1962). 
[13] L. Belloni, Phys. Rev. Letts 57, 2026 (1986). 

[14] A. D. Kirschembaum, J. A. Cahil, P. J. McGoingle, and A. V. Grosse, J. Inorg. Nuclear Chem. 
24, 1287 (1962). 

[15] P. N. Vorontsov-Vel'yaminov, A. M. El'yasevich, L. A.Morgenshtern, and V. P. Chasovskikh, 
Teplofiz. Vys. Temp. 8, 277 (1970). 

24 



[16] P. N. Vorontsov-Vel'yaminov and V. P. Chasovskikh, Teplofiz. Vys. Temp. 13, 1153 (1975). 
[17] J. M. Caillol and J. J. Weis, J. Chem. Phys. 102, 7610 (1995). 

[18] E. Luijten, M. E. Fisher, and A. Z. Panagiotopoulos, Phys. Rev. Lett 88, 185701 (2002). 
[19] A. Z. Panagiotopoulos, J. Chem. Phys. 116, 3007 (2002). 

[20] J. M. Romero-Enrique, G. Orkoulas, A. Z. Panagiotopoulos, and M. E. Fisher, Phys. Rev. 

Lett 85, 4558 (2000). 
[21] Q. Yan and J. J. de Pablo, J. Chem. Phys. 114, 1727 (2001). 
[22] Q. Yan and J. J. de Pablo, Phys. Rev. Lett. 86, 2054 (2001). 
[23] A. Z. Panagiotopoulos, Phys. Rev. Lett 88, 45701 (2002). 

[24] J. M. Romero-Enrique, L. F. Rull, and A. Z. Panagiotopoulos, Phys. Rev. E 66, 041204 (2002). 

[25] Q. Yan and J. J. de Pablo, J. Chem. Phys. 116, 2967 (2002). 

[26] H. L. Friedman and H. Taube, J. Am. Chem. Soc. 72, 3362 (1950). 

[27] M. L. Japas and J. M. H. L. Sengers, J. Phys. Chem. 94, 5361 (1990). 

[28] R. R. Singh and K. Pitzer, J. Chem. Phys. 92, 6775 (1990). 

[29] H. Weingartner, S. Wiegand, and W. Schroer, J. Chem. Phys. 96, 848 (1992). 

[30] T. Narayanan and K. Pitzer, J. Phys. Chem. 98, 9170 (1994). 

[31] E. Gonzalez- Tovar, M. Lozada-Cassou, L. M. y Teran, and M. Medina-Noyola, J. Chem. Phys. 

95, 6784 (1991). 
[32] M. Lozada-Cassou, J. Chem. Phys. 75, 1412 (1981). 
[33] M. Lozada-Cassou, J. Chem. Phys. 80, 3344 (1984). 
[34] M. Lozada-Cassou and E. Di'az-Herrera, J. Chem. Phys. 92, 1194 (1990). 
[35] M. Lozada-Cassou and E. Di'az-Herrera, J. Chem. Phys. 93, 1386 (1990). 
[36] J. E. Sanchez and M. Lozada-Cassou, Chem. Phys.Lett. 190, 202 (1992). 
[37] M. Lozada-Cassou, in Fundamentals of inhomogeneous fluds, edited by D. Henderson (Marcel 

Dekker, New York, 1993), Chap. 8. 
[38] J. Alejandre, M. Lozada-Cassou, E. Gonzalez- Tovar, and G. A. Chapela, Chem. Phys. Lett. 

172, 111 (1990). 

[39] J. P. Valleau, R. Ivkov, and J. M. Torrie, J. Chem. Phys. 95, 520 (1991). 
[40] M. Lozada-Cassou, W. Olivares, and B. Sulbaran, Phys. Rev. E 53, 522 (1996). 
[41] All those theories which are derived from the Orstein-Zernike equation plus any of the estab- 
lished closures, we refer to as standard integral equations theories. 

25 



[42] R. Kjellander and S. Marcelia, Mol. Phys. 83, 789 (1994). 

[43] M. Plischke and D. Henderson, J. Phys. Chem. 92, 7177 (1988). 

[44] M. Plischke and D. Henderson, Electrochim. Acta 34, 1863 (1989). 

[45] P. Attard, J. Chem. Phys. 91, 3072 (1989). 

[46] P. Attard, J. Chem. Phys. 91, 3083 (1989). 

[47] M. E. Fisher, J. Stat. Phys. 75, 1 (1994). 

[48] Y. Levin and M. E. Fisher, Physica A 225, 164 (1996). 

[49] N. Bjerrum, Kgl. Dan. Vidensk. Selsk. Mat.-Fys. Medd. 7, 1 (1926). 

[50] B. Hafskjold and G. Stell, in The liquid state of matter, edited by E. W. Montroll and J. L. 

Lebowitz (North Holland, New York, 1972). 

[51] G. Stell, K. C. Wu, and B. Larsen, Phys Rev. Lett 37, 1369 (1976). 

[52] Y. Q. Zhou, S. Yeh, and G. Stell, J. Chem. Phys. 102, 5785 (1995). 

[53] S. Yeh, Y. Q. Zhou, and G. Stell, J. Phys. Chem. 100, 1415 (1996). 

[54] F. O. Raineri, J. P. Routh, and G. Stell, J. Phys. IV 10, Pr5 (2000). 

[55] L. Blum and O. Bernard, J. Stat. Phys. 79, 569 (1995). 

[56] O. Bernard and L. Blum, J. Chem. Phys. 104, 4746 (1996). 

[57] J. W. Jiang, L. Blum, and O. Bernard, Mol. Phys. 99, 1765 (2001). 

[58] J. W. Jiang et al, J. Chem. Phys. 116, 7977 (2002). 

[59] S. H. Suh, L. M. y Teran, and H. S. W. H. T. Davis, J. Chem. Phys. 142, 203 (1990). 

[60] K. C. Zhang, M. E. Briggs, R. W. Gammon, and J. M. H. L. Sengers, J. Chem. Phys. 97, 
8692 (1992). 

[61] S. Wiegand et al, J. Chem. Phys. 109, 9038 (1998). 

[62] R. R. Singh and K. S. Pitzer, J. Am. Chem. Soc. 110, 8723 (1988). 

[63] T. Narayanan and K. S. Pitzer, Phys. Rev. Lett. 73, 3002 (1994). 

[64] T. Narayanan and K. S. Pitzer, J. Chem. Phys. 102, 8118 (1995). 

[65] K. Gutkowski, M. A. Anisimov, and J. V. Sengers, J. Chem. Phys. 114, 3133 (2001). 

[66] D. A. McQuarrie, Statistical Mechanics (Harper and Row, New York, 1976). 

[67] H. L. Friedman, A Course in Statistical Mechanics (Prentice Hall, Englewood Cliffs, NJ, 1985). 

[68] J. P. Hansen and I. R. McDonald, Theory of simple liquids, 2nd ed. (Academic Press, London, 
1986). 

[69] M. Born and M. S. Green, A General Kinetic Theory of Liquids (Cambridge University Press, 

26 



Cambridge, 1949). 

[70] J. Yvon, La Theorie Statistique des Fluides et I'Equation d'Etat, Actualites Scientifiques et 

Industrielles (Hermann, Paris, 1935), Vol. 203. 
[71] G. Arfken, Mathematical Methods for Physicists (Academic Press Inc., San Diego, CA, 1985). 
[72] R. Messina, C. Holm, and K. Kremer, Phys. Rev. Lett. 85, 872 (2000). 
[73] R. Messina, C. Holm, and K. Kremer, Europhys. Lett. 51, 461 (2000). 
[74] R. Messina, C. Holm, and K. Kremer, Phys. Rev. E 64, 021405 (2001). 
[75] M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 (1998). 

[76] L. Degreve, M. Lozada-Cassou, E. Sanchez, and E. Gonzalez- Tovar, J. Chem. Phys. 98, 8905 
(1993). 

[77] In Fig. 4 of ref. |2J it is plotted the energy calculated by the unscreened interaction potential 
of four ions, and is not a simulation result 8(|. In our case, our conclusions are derived from 
the potential of mean force. Hence aligned clusters seem to be better supported by our theory 
and is not in contradiction with the simulation results of Yan and de Pablo. 

[78] L. M. y Teeran, E. Diaz-Herrera, M. Lozada-Cassou, and R. Saavedra-Barrera, J. Comp. Phys. 
84, 326 (1989). 

[79] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, 4th ed. (McGraw Hill Book 

Company, London, 1989), Vol. 1. 
[80] J. J. de Pablo, personal comunication. 



27 



n Exact 


HNC 


PY 


■ A 


A 


A 


2 mm 







FIG. 1: Mayer diagrams for the first (n = 1) and second (n = 2) order in p expansion of the 
pair correlation function, g^ 2 \r). The exact, hypernetted chain (HNC) and Percus-Yevick (PY) 
coefficients are shown. The black points and white dots are called field and root points, respectively. 
The bonds represent the Mayer function /(ru)- 

particle 3 




FIG. 2: Schematic representation of the model. 
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FIG. 3: Two examples of the transformation of Mayer diagrams under TPE. The notation for the 
particles and species is the same as in Fig. |2I The dashed bond represents the f 1 i{r2z) = d "j r ^ 23 ^ 
function. Thereby, in (N,5), N stands for the particle number and 5 for the particle species (see 
also Fig. |2J). N = stands for the dumbbell particle, (a) An example of a second order Mayer 
diagram involving a dumbbell particle and particle 3. (b) The same diagram as in (a) but the 
constitutive particles of the dumbbell are explicited, i. e., at the level of the triplet correlation 
function, (c) Resulting diagrams upon applying BGY (on the level of the mean force). 
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FIG. 4: Three particle distribution function g ++i (r,6 = 7r/2;r = a) for a dumbbell made of two 
(divalent) positive particles, with p = 1M and z = 2. The solid lines represent the results from 

TPE-HNC/MSA. The MD results are shown in filled and open circles for g+\ (r,ir/2;a) and 

(?+++(?", vr/2; a) respectively. 
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p = 1M ; ++ ; 6 = 71/4 




FIG. 5: Same as in FiglUwith 9 = 7r/4. 
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FIG. 6: 3D representation (in Cartesian coordinates) of: (a) (upper 3D plot) #++_(r, #;r = a) and 

[31 

(b) (lower 3D plot) g^ + _(r,6;T = 2a) obtained by TPE-HNC/MSA for the same fluid parameters 
as in Figs. |1] and |SJ The dumbbell axis is parallel to y axis. 
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FIG. 7: Three particle distribution function g + _i(r, 9 = 7r/2;r = a) for a dumbbell made of a 
positive and a negative divalent ions, with p = 1M and z = 2. The solid lines represent the results 
from TPE-HNC/MSA. The MD results are shown in filled circles. The curves for g^L_(r, 7r/2; a) 

[31 

and g + _ + (r, 7r/2; a) colapse in a single curve. 
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p = 1M ; +- ; 6 = 71/4 
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FIG. 8: Same as in Fig. E| with 9 = vr/4. The MD results are shown in filled and open circles for 

f3l [3l 

g + (r, 7r/4;o) and g + _ + (r, 7r/4; a) respectively. 
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p= 1M;2:2 

10 i — . — . — . — . — i — . — . — . — . — i — ^ 




-10 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 
1 2 3 4 5 

r/a 



FIG. 10: Effective forces between two like charges [F ++ (r) > 0] and two opposite charges 
[F^ (r) < 0] as a function of their separation r, in reduced units of — — , with p = 1M and 

IB 

z = 2. The solid lines represent the results from TPE-HNC/MSA and the dashed lines are the re- 
sults from HNC/MSA. The MD results are shown in filled and open circles for (r) and F ++ (r), 

respectively. 
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FIG. 11: Three particle distribution function with p = 0.005M and z = 2: (a) g+n i r i 7r /2;a) and 

(b) 5 [ f 3l f+ (r,vr/2;a). The solid lines represent the results from TPE-HNC /MSA. The MD results 

are shown in filled and open circles for (a) g^} (r,ir/2;a) and (b) g^_ + (r,ir/2;a), respectively. 
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FIG. 12: Same as in Fig. ^2 with # = 7r/4. The MD results are shown in filled and open circles 
for g+\ (r,n/2;a) and <jf!j!!|_ + (r, 7r/2; a), respectively. 
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FIG. 13: 3D representation (in Cartesian coordinates) of g + \ (r, 9; r) obtained by TPE-HNC/MSA 

for the same fluid parameters as in Figs. Illl and fT^l The dumbbell axis is parallel to y axis, (a) 
r = a (b) r = 5a. 
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FIG. 14: Three particle distribution function g + _i{r, 7r/4; a) with p = 0.005M and z = 2. The solid 
lines represent the results from TPE-HNC/MSA. The MD results are shown in filled and open 
circles for g + (r, 7r/4;a) and g + _ + (r, it/A; a), respectively. 
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FIG. 15: 3D representation (in Cartesian coordinates) of g , (r, 9; a) obtained by TPE-HNC/MSA 

for the same fluid parameters as in Fig. 1141 The dumbbell axis is parallel to y axis. 
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FIG. 16: Effective forces between two like charges [F ++ (r) > 0] and two opposite charges [F-\ (r) < 

0] as a function of their separation r, in reduced units of — - — , with p = 0.005M and z = 2. 
The solid lines represent the results from TPE-HNC/MSA and the dashed lines are the results 

from HNC/MSA. The MD results are shown in filled and open circles for (r) and F ++ (r) , 

respectively. 
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FIG. 17: An example of the grid (in Cartesian coordinates) used to solve Eq. l)Aljl . 



N Finite Elements 




FIG. 18: The same grid as in in Fig. El but mapped into the rj — £ plane. A triangular element 
and its 6 nodes are represented. 
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